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Abstract 

In this paper, we apply Value-At-Risk (VaR) approaches on the problem of yearly electric 
generation management. In a classical approach, the future is modelled as a markov chain and 
the goal is to minimize the average generation cost over this uncertain future. However, such a 
strategy could lead to big financial losses if worst case scenarios occur. The two VaR approaches 
we propose, precisely aim at robustifying the model. On a practical point of view, it amounts to 
introduce a new set of constraints modelling the uncertainties in the original optimization problem 
or equivalently to change the dual objective function. The new optimization problems are solved 
as efficiently as the nominal model. Numerical simulations are presented and discussed for this 
application. 

1 Introduction 

In this paper, we are interested in optimization problems arising in electrical yearly power 
management. Given electric generation plants (nuclear, thermal and hydroelectric power generator 
plants, demand side management contracts modelled as a virtual plant called E.JP), the objective is 
to minimize the production cost over a yearly horizon, to fulfill operating constraints of generation 
units and the equilibrium between production and demand at each time step. In practice, the 
modelling approach is highly depending on the time horizon of the optimization problem : for 
short time horizons, typically daily or weekly, the problem is generally assumed to be deterministic 
(cf. n,(2l), whether for longer management horizons, a special emphasis is done on the stochastic 
nature of data and events. In particular, on a yearly scale, reservoir inflows, demand, availability 
of the plants as well as electricity prices cannot be considered to be deterministic : in France, 
for example, winter costumer's demand has uncertainties that can reach one GW per decreasing 
temperature degree while the peak loads are around 70 GW ! So a true challenge is to ensure 
robustness of computed optimal production and marginal value face to various uncertainties like 
customer demand but also water inflows or plants unavailability. 

In general, for yearly generation management, utilities are not interested in generation schedul- 
ing but rather management strategies or Bellman values to perform Monte Carlo analysis of the 
futur. Since Stochastic Dynamic Programming quickly comes to its limits for the optimization of 
high dimensional state systems, a decomposition approach is usually necessary [H]. In our case, 
a large scale numerical optimization is first formulated and solved solved using Lagrangian relax- 
ation to provide marginal costs on a scenario tree. Then, those marginal costs are used to compute 
local feed-backs (see below section 2.1 for details). This decompostion framework that could be 
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compared to the Dual Stochastic Dynamic Programming method Q is based on the adaptation of 
[H] to yearly generation management. The main drawback of this scheme relies in the local aspect 
of the feed-back functions that loses the robustness of the global Bellman function. 

In this paper, we show that a Value- At-Risk (VaR) approach for modelling uncertainties allows 
to enhance the robustness properties of those local feed-back laws with the following practical 
benefits: (i) significant reduction of the variance of simulated cost - up to 38% of reduction for a 
comparable average cost if the worst case scenario occurs -; (iii) parsimonious use of water reservoir; 
(ii) reduction of very high cost strategies. Moreover, we will show that there exits VaR approaches 
on the dual optimization problem that preserve the space decomposition approach while having 
a very nice economical and physical interpretation. The first one comes from a primal relaxation 
on the demand side - how to control the sales turnover face to uncertainty of the demand ?- and 
the second one from a dual relaxation -how to control the production costs face to plant random 
unavailability ?-. 

The sequel is organized as follows. In the next section, the optimization model is described and 
special focus is given on random events and inputs. Section 3 deals with the robustification issues, 
the VaR approaches, the connection between VaR and duality as well as implementation issues. 
Section 4 presents some numerical results comparing the performances of the nominal and robust 
models and finally some additional modeUing details are given in the appendix. 

2 Setting of the physical model 

The physical model is a stochastic dynamic system where the random inputs are the costumer's 
demand, the unavailability of the thermal units and the quantity of natural water infiows. With 
a representation of the random inputs and events as Markov chains, we can naturally formulate 
the optimization problem as a finite horizon discrete time stochastic control problem on a scenario 
tree representing the behavior of the random inputs and states. 

2.1 Model setting 

We aim at minimizing the average production cost along the scenario tree. If we describe 
at each node n of the tree the states of the plant production unit £ by the variable and the 
commands applied to this plant by u^, the stochastic control problem may be formulated as ||6J: 

eec neo 

where 

• £ is the set of plants and O is the set of the nodes, 

• 7r„ is the probability to be at node n, 

• Vn is the set of time subdivisions associated to node n, 

• {p) the state of plant £ at node n and time subdivision {p) , 

• ul^{p) is the control variable of plant i at node n and time subdivision (p), 

• Pn,pi^nip)>''^n,pip)) the production of plant £ in the state xf^ when command u^{p) is 
applied to this plant at node n and time subdivision (p), 

• Vnip) is the costumer demand at node n and time subdivision (p). 
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• is the production cost when command u^{p) is appUed to unit £ in the 
state xl^{p), 

• Xe is the functional set of constraints on the control and state variables of plant £. 

In fact, due to the autonomy of the plants, the model may be reformulated as a linear op- 
timization problem with separated domains of constraints and one coupling constraint (produc- 
tion/demand equilibrium). Each domain of constraint is a dynamic system describing a plant 
process. Introducing cost vectors Ci and coupling matrices Ai with ad hoc sizes for i £ {1, ...,4}, 
we can formally describe Q in the following way: 

min c[ut + c^Un + cT^Xh + cjxe 

{ut,Un, Uh, Ue) € T X Af X H{xh) X f (Xg), (2) 

AlUt + A2Un + A^Uh + AaUe ~ d e M-°, 

where each control variable u, belongs to a set parameterized^ by the corresponding state variable 
X, : 

• ut, control variable of classical thermal plants, with T as the set of constraints for the thermal 
plants subset; 

• Un, control variable of nuclear thermal plants, with J\f as the set of constraints for the thermal 
plants subset; 

• Xh , state variable of hydrauHc plants with Uh as the control variable and dynamics described 
by Uh e H{xh); 

• Xe, state variable of EJP contract, with Ue as the control variable and dynamics described 
by Ue e Sixe); 

• d £ , the vector of demands with D = (n° of time subdivision ) x card (O). 
Therefore we may write ^ as: 

{min f^(u, Xh, Xp) 
ueu^ ' (3) 
Au = d given e R-^ , 

when setting 

U = {ut,Un,Uh,Ue) e U = Tx J\f X H{xh) X £ (Xe) , (4) 
f{u, Xh,Xe) = C^Ut + C^Un + C^Xh + C^Xe, (5) 
Au — AlUt + A2Un + A'iUh + A/^Ue. (6) 

Notice that in this model, the demand c? is a fixed vector corresponding to the reaHzations of the 
demand at the different nodes to the scenario tree. 

2.2 Model Analysis 

The efficiency of this model is assessed on a set of independent scenarios representing different 
evolutions of the demand, the inflows for hydro reservoirs and the outages of the thermal units. For 
each scenario, a generation schedule as well as its cost are determined (a detailed description of the 
implementation of the generation schedule is given in the appendix) . Such a model has intrinsic 
limitations essentially linked to the fact that the strategy computed by the algorithm above is 
optimal only on the optimal trajectory of each reserve. In this sense, the Bellman functions 
obtained by Stochastic Dynamic Programming on marginal values only give a local optimum. The 

^see Appendix A2 for a detailed plants description 
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effect of using such local feedbacks as global strategies is to create a high volatility of scenarios 
costs. So it is desirable to strengthen this model by including a more reliable model of uncertainty 
on the scenarios at the earliest stages of the problem setting. The goal is to ensure some regularity 
of the optimal strategies with respect to the inputs of the optimization problem. In other words 
we would like to find a robust counterpart with the following properties: 

(i) reduce the volatility of the simulated costs over a continuum set of reasonable scenarios; 

(ii) reduce the number of extreme case optimal strategies (parsimonious use of water reservoir 
that might not be nearly empty for a long period); 

(ii) reduce the number of very high cost optimal strategies. 

We will see in the next section that such an objective of variance reduction may be easy to formulate 
in a Value- At-Risk setting. 



3 Robust Counterpart of the Decision Model 

3.1 The Value- At-Risk approach for stochastic optimization problems 

Let {(jj,r,x) — *■ f{r{uj),x) be a concave (with respect to x) income functional depending on a 
random function lu r{uj) where a; £ X C M" is deterministic variable, X being a non empty 
closed and bounded set. A Value- At-Risk (VaR) approach allows us to choose x leading to the 
maximal possible income with a given confidence level < e < 1. Typically, if we have additional 
constraints on x expressed as g{x) > 0, one formulates the following optimization problem : 

max 7 

P{f{r{w),x)>j)>l-e, (7) 
X £ X, g{x) > 0. 

Let $ be the cumulative distribution function of the Gaussian density. Following [3 and [H] , we 
can find a Risk Averse solution (indeed a prudent one as upper bound of the optimum) in some 
usual cases when solving: 

CT/^p^/"^^^ ^Miri^^),x)]- K{e)a^[f{r{u;),x)] , 
\x£X, g{x) > 0, 

where K{e) is a risk factor depending on the assumptions on the distribution: 



$ ^(1 - e) > if /(r(»), x) is Gaussian, 
^if/(r(.),x)eLi nLi 



<^) = { fT^ -.c ^^ ^ rl r.r2 (9) 



For instance is f{c{ijj),x) — c{lo)'^x is a linear function of x, the problem {VaRe) simply reduces 
to 

max E^[c{ujY']x — K{e)V x'^Tx with — cov{ci{w), Cj{w)), 
x£X, g{x) > 0, ^ ' 

and if F is invertible the above problem appears as the Robust Counterpart of the problem 

max c{uj)'^x . ^ 

x£X, g{x) > 0, ^ ' 

where the uncertainty set chosen for the random vector c{uj) is the ellipsoid: 

{x £ E", {x - E^[c{iu)]fT-\x - E4c{lu)]) < K^{e)}. 
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Now it is clear that a VaR approach is a practical way to calibrate a variance penalty term for 
a maximization of a random functional. In the case of power generation management, it means 
that we aim to find the best compromise between production cost and volatility of strategies at 
(possibly) the extra cost of some sub-optimality on the most favorable scenarios. This approach is 
very easy to set as a regularization approach for linear programs because one only needs to specify 
a risk exposure level e. 

Remark 3.1 Being a covariance matrix, F is symmetric and positive semidefinite, and using the 
induced norm, the objective function in ifTIill may be written as E^[c{uj)x] — K(e)||a;||r. Therefore 
a VaR regularization appears as an Han Penalization for the problem: 

max E^[c{u;,x)] 

xeX, g{x)>0, (t[c{uj,x)] < a, ^ ' 

where the penalization coefficient is chosen on a probabilistic basis. We mention that in the convex 
case, there exists some sufficient conditions linking K{e) to the dual norm of || • ||r that ensures 
the exactness of this penalization. Remarking that if F ;^ 0, the dual norm is || • ||r-i; we can 
mention that the condition given in |9' reduces to K{e) > ||A||r-i where A is the lagrange multiplier 
associated to the optimal solution. Therefore : 

(i) a sufficient condition for exact penalization is e < Yq:x~(ry' where Ami„(r) is the smallest 
eigenvalue of F. 

(ii) Optimizing the smallest eigenvalue of the covariance matrix F will allow to reduce the bound 
on £ and then to enforce the constraint of risk reduction: P{c{uj)x > j) > I — e. So any 
matrix calibration technique that will reduce the condition number of the covariance matrix 
by increasing the smallest eigenvalue will provide additive degree of freedom if we use a 
probabilistic constraint to control the income. 

3.2 Application to the Power Generation Model 

The idea is to take advantage of a decomposition of the dual optimization problem to introduce 
a VaR modelization on two subproblems : first on the uncertainty on the demand and next on the 
unavailability of thermoelectric plants. To begin with, we point out that the dual problem of lO 
is max 0(A) where: 

e{\) = 0d{X) + 0t{\) + 0n{\) + 0HiX) + 0j{X) - 0d{X) + 0{X), (13) 

once we have introduced the partial dual functions: 

0d{X) - X^d, 
OriX) - inf (ci - A^Xf ut, 

0n{X) = Jl^^ {C2 - A2 X)'^ Un, ^^^^ 

0h{X) = inf c^Xh - X^AsUh, 

0j{X) — inf C^Xe — X^AiUe- 

lie ^£{Xe) 

Now it is natural to robustify ijSJ by formulating a VaR problem on the subsystem with random 
information or state : 9d for the costumer's demand d and {0t, 0n,0h) for the states of the thermal, 
nuclear and hydro plants. 

3.2.1 Primal Relaxation and Dual VaR regularization of the demand 

The idea is to make a primal relaxation on the predictions made for the demands on the different 
scenarios (i.e the values of the demand d at the different nodes of the scenario tree) that are prone 
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to errors. Rather than considering that the demands at each node and each time subdivision of 
the tree are known exactly, we suppose that d belongs to a given uncertainty set £ which is the 
ellipsoid given by: 

£ = £{d,r,K) = {xGM.^ I (x-dfY-^(x-d)<K^{e)}. (15) 

where d ^ Ej{d{w)\ the covariance matrix F is given by = cov{di{oj), djUj)) and K{e) depends 
on the assumptions made on the distribution of the demand (see HH]). That means that we 
reformulate the problem JSJ as: 

{min f^(u,Xh,Xe) 
Au = d e £. 

This is a relaxation of problem (EJ: the demand vector is no longer fixed but can be any vector 
from the eUipsoid £. Solving ltT6|l by duality amounts to solve max07j(A) with 

^TzW — niin f^{u, Xh, Xe) + {d — Au) 

= eT{X)+ON{X)+OH{X)+Oj{X)+mmX'^d. ^ ' 

Now notice that min A"^d = 0£(A) where is the support function of the uncertainty set £ and is 

deS 

given by : 

ct)£{\) ^ X^d- K{e)VX'^T\. (18) 

Note that the robustification just turns out to replace 9d{X) = X^d in the original dual function 9 
by (t>£{X) — X^d— kVVFA, where d = Eij[d{uj)]. We can notice that the relaxation of the demand 
in the ellispoid £ in problem lll(i|l amounts to use a VaR approach on the dual problem of problem 
(LP). Indeed, using a VaR orientated technique, as the demand d is random, instead of maximizing 
0{X) (which is the dual problem of problem (LP)) we could maximize 0{X) +7*(A), where 

f max 7 

^ (^) = | P(AM-)>7)>l-e. (^^) 

From subsection Kill this VaR approach reduces to max 6ti{X), which is the dual problem of prob- 
lem {LP-ji). In what follows, this VaR approach will be denoted by VaRpA- 

Economical Interpretation. The problem Ijl9|l may be interpreted as the maximization of 
the minimal Sales Turnover that can be ensured with an arbritary degree of confidence. In other 
words, a performing regularization strategy by a relaxation of the costumer's demand leads to the 
reduction of the volatility of the sales turnover. 



3.2.2 VaR approach on the dual thermal problem 

In this subsection, we intend to exploit a stochastic model of the unavailability of the ther- 
mal plants in order to formulate a Value- At-Risk problem on the costs of thermoelectric power 
generation. Let £ be a thermal unit with thermal groups. Let aj^e{t) be the probability that 
group j of unit £ works at time step t and C/* ^ the random variable such that — 1 if group j 
works at time step t and Uj ^ = else. We suppose that the groups are regularly checked and, if 
necessary, repaired every toq time steps. Between two consecutive checking dates, we assume that 
the availability of the units is not changing. This means that between two consecutive checking 
dates, a given group is either working or it is out of work during the whole period. If io = 1 and 
tk = mok for k G N* , then the probabilities aj,i{t), for t ~ tk, ■ ■ ■ , tk+i — 1 are the same and we 
only need to evaluate aj^i{tk), k >0. Those probabilities aj^i{tk) that a group j of unit £ works at 
time step mok will depend on the past evolution of the availability of this group. If at time step 
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tk-1, the group was out of work, there is a big probabiHty (say I — (3f with /3f small) that it works 
at time step tk (the time between two checking dates is greater than the mean time to repair) and 
a small probability /?f that it is still out of work at time step tk . Now if the group was working 
for the last m periods delimited by the last m + 1 checking dates, we can assume that the longer 
it has been working without failure (the larger m) the more likely it can break down at time step 
tk- Thus, there is a decreasing function of m, /?2("i) such that for any group j of unit £, 

P{U'j} = 1 1 Group j was working from tk~m to ife-i) = (32{m). 

A particular case is the case where the state process of a given group is an homogeneous Markov 
chain where the state space is {F,W} where F stands for the failure state and W for the working 
state. In this case, (32{m) = /?| is fixed and corresponds to the probability for a group of unit £ to 
work on a given period knowing that it was working the period before. The transition matrix for 
the groups of unit £ is given by: 



1-/3^ 

The probability aj^i{tk) is then given for fc > 1 by: 

= Pf(j)^^/(1, 2) + pV(j)^/(2, 2); 

where p^{j) = 1 — Ppij) — and p^{j) is the probability that group j of unit £ works at the 
first time step. For the simplicity of the exposure, we assume that for a given unit £, either all 
the groups are working or all the groups are out of work at the first time step. Thus, aj^e{tk) 
is j-independent and ag{tk) will denote the probability that a group of unit £ works at time tk- 
Further, we can partition the scenario tree in subtrees such that the root node and the leaves nodes 
of a given subtree respectively correspond to time steps tk and tk+i for some fc G N. Thus the 
unavailability rates at the different nodes of any subtree of this partition are the same for a given 
unit. Let O = iik=i Ok be such that Ok are the nodes of the k-th subtree Sk in this partition. Let 
Tk = {{jiP), I j & Ok,p ^Vj}- The dual thermal subproblem then writes: 

/ in 

min {cuk ~ >^kY' uuk 

I k=l 

< Utik < n{k)T^{k)P^^^dk, 

where Pj^ax is the maximal available power of thermal unit £, uuk = i'^tejp)(j.p) e t^j '^uk = 
{cuip){j,p) GTfc, Afe = {Xj,p){j,p) (ZT^, dk = {dj,p){j.^p) 6Tfc (see the appendix), T-f(fc) gives the pro- 
grammed unavailability rates for unit £ and the times subdivision of the set Tk and r^(fc) is the 
unavailability rate of unit £ for the nodes of the set Ok- We reformulate this problem as the 
following problem: 

C m 

min ^^T£(fc) {ciik - Xk^utik 
e k=i 

< Utik < T^{k)P^^^dk, 

by setting uuk = ^ffk)- Notice that for quite a number of linear stochastic optimization problems, 
the random is only in the right hand side of the constraints as it is the case for the dual thermal 
subproblem. The above simple transformation allows to transfer the random in the objective 
and to implement a VaR method to compute robust solutions as was described in subsection l3.1l 
Given a confidence level < e < 1, we thus now introduce a VaR approach on the thermal plant 
cost/revenue balance: 

min 7 

{VaR)Benef l ^EH^^^^) (^"^ ~ ^kf ^Uk <l)>l-S, (20) 

e k 

< uttk < T^{k)P^^^dk- 
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This problem may be understood as a problem of maximization of the benefits or equivalently 
as a problem of minimization of the losses. We now need to study the modelling of the un- 
availability rates Ti{k) to give an explicit form for problem Ij2()|l . 

Modelling of the unavailability rates Ti{k). Let P,^ax be the maximal power of a group in unit 
£. Then the theoretical maximal power available on the thermal unit t is given by -P^ax = '^^^max- 
The maximal power available of unit t for the nodes of the set Ok is then: 



rii y^«f 7-7-tfc 

-'max ~ / , '-^ j.t mia^ ~ max „ ~ max ' ^ l,"-^ ■ 



^ - m 

If t{k) is the time step associated with the root of the subtree Sk, notice that under the above 
hypothesis, the random variable npTii{k) follows the binomial law B{n(^,ai{t{k))). We then have 
E[n{k)] = ai{t{k)) and the variance oinik), \&x(n{k)) = "^(«(fc))(i^-"^(*(fc))) _ Xi^k be the 

random variable reXk) (cuk — ^k)'^utik- As the {Ti{k))i^k are independent we have: 

i,k 

E[Xt^k\ ^ oieXt{k)){cuk - \k)^utik, and VaR[Xi^k\ = ^ Qtk utik, 
e,k e,k 

where the matrix Qik is defined by (^cigk — ^k){cuk ~ Afe)"^. From subsection 13. 11 

ll2?11l amounts to solve: 



min y^ai(t{k)) (cuk- >^k)^uuk + n(£) ly^uJgkQikUuk 

{VaR)Benef{ it ^7% (21) 

< uuk < 4{k)Pi,,dk, 

where K{e) — \f^^- Notice that this is a second order cone optimization problem whose dimension 
will be high in practice (the number of nodes of the tree times the number of subdivision times). 
Another more conservative approach would be to use a VaR approach for each time subdivision of 
the nodes. This is possible as the dual thermal subproblem is separable with respect to the time 
subdivisions of the nodes. The advantage of this approach is that we have an expHcit solution for 
the new dual thermal subproblem and the new dual thermal subproblem is again separable with 
respect to both the thermal units and the time subdivisions. This allows to solve problems of big 
sizes. Moreover, the unavailability rate of unit £ for each time subdivision follows a binomial law 
which can be approximated by a Gaussian law if {ng ai{t{k)) > 10 and (1 — ai{t{k))) > 10) or ni 
big enough, say ng > 6. At last, we mention that a particular case would be to consider that ae{tk) 
doesn't depend on k. In what follows we both suppose that ai{tk) doesn't depend on k and that 
the VaR approach on the thermal subproblem is done at each time subdivision. This will allow to 
check how the VaR approach on the thermal subproblems gives an immunization with respect to the 
uncertainty we have on the unavailability of the units which works surprisingly very well in practise. 

Economical Interpretation. Qualitatively, one can say that the objective of (VaRBenef) is 
to maximize the benefits of a thermal unit while ensuring some robustness with respect to plants 
random unavailability. 

3.2.3 Intermediate summary 

At this point, we have introduced two different regularizations on the original dual optimization 
problem that was initially formulated as : 

(DP) max 9d{X) + max {Ot + On + Oh + dj) (A) (22) 
where the 6j are given by dJ. 
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1. From a relaxation of the demand, we formulate the dual regularized problem 

{VaRpA) msx9f{\) + max (6It + 6Ijv + 0h + Oj) (A) (23) 

with 

ef{X) = E^{d{uj)Y\-K{e^)VWvX. (24) 
where ei is the confidence level chosen to implement VaRpA- 

2. From a rewriting of the random unavailability of the thermal plants, we get 

[VaRBenef) max6ld(A) + max6'P-(A) + max(0Ar + 6*// + Oj) (A) (25) 



with 



n,p,£ECT ■ 

Vl(u, Xn,p) = (^ae{cunp - Xn,p) + K{£2)\f^^^^^\cUnp - An,p|^ 



(26) 



where 62 is the confidence level chosen to implement VaRBenef ■ If the we use a VaR technique on 
the whole thermal subproblem then 9^{X) is given by ll2T|l . Finally, we can also combine the two 
previous regularizations as the mixt problem: 

(Fai?™„t) max e'^-JX) = max0j(A) + max0?(A) + max(0jv + 0h + 9.,) (A). (27) 



4 Implementation and Numerical simulations 

4.1 Implementation and simulation protocol 

To solve the primal optimization problem ijSJ, we have to solve minmaxL(u, Xh, Xe, A), where L 

is the usual Lagrangian. This will be equivalent to the dual problem max 9{X) if only thermal and 

hydro units are considered. Indeed, in this case, problem ijSl is a below bounded linear program 
and both the primal and the dual are equivalent to each other. If we take into account EJP 
contracts, the set of constraints is not convex and the duality gap is strictly positive. However, 
the weak duality relationship still holds: 

minmax L(u, a;^, Xe, A) > max 6{X). (28) 

u^ti A A 

Moreover, numerical simulations have shown that the duality gap is generally quite small. The 
dual problem thus allows us to approximate primal solutions and estimate marginal prices. 



4.1.1 Space decomposition for Optimization 

First, we describe the space decomposition method for the dual function 9 and next, we explain 
the adaptations necessary for the regularized problems. The dual function 9 is non differentiable, 
concave and separable with respect to the units as it writes 9(X) = X^d + ^ ^^(A) with 

eec 

= , min V V 7r„C^ (x^(p),ul(p))-A„,pP^ (a;l(p),u^(p)), 

the dual function of the subproblem associated with unit I. This is especially of interest to treat 
problems of big size as it is the case for our application. To maximize 9 (or which is the same to 
minimize the convex function —6') we use a bundle method described in ^l]. This requires to build 
a black box which, for any A e MP is able to compute — ^(A) and to give an arbitrary subgradient 
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s(A) e d{—6{X)). Let Lg be the partial Lagrangian associated with the partial dual function 
corresponding to unit £. The computation of —0{X) is done solving the different optimization 
problems associated with the different production units. As for a computation of a subgradient, if 
/(A) = argmin^_„ A) and P^(A) = {Plpivi^pW))^^^ then 

-d\fi) > -0\\) + {pL- Xf P\\) for all (A, /^) G R^, 
which shows that -P^(A) G d{—9^{X)). We then immediately have 

siX) = -d + Y, P'iX)edi-9iX)). 

More precisely, following ?9l, a global resolution by an iterative scheme can be described with 4 
steps starting from a reference price A used to initialize the algorithm with index k = 1 and Ai — X. 

1. At iteration k, decomposition in subproblems and computation of the local solution of sub- 
problem £: y^(Afc); 

2. Evaluation of the dual function 6 at the point Afc and computation of a subgradient s(Afc); 

3. Updating of the multipliers by the coordinator using a black box method (i.e computation 

of Afc+i); 

4. Updating of the index k ^ k + 1 and go to step 1. 

The robustifications proposed still remain within the same framework, so we can solve them in 
a space decomposition framework. The adaptation of the above algorithm to the robustifications 
proposed is easy since we just need to modify step 2. 

• For the dual regularized problem {VuRpa), we have to increase the value of the dual function 

YX 

by —K{ei)V X^rX and the sugradient by — k(£i) ^tyx ^'^^^^ confidence level used 

to implement {VaRpA)- So the only extra cost of this model is the estimation of matrix F. 

• The method VaRbenef simply modifies the thermal dual problem which becomes problem 
(j2()ll (if the VaR approach is done on the whole tree) or Ij2(ill else. If (j2(ill is used, then the 
solution of the new thermal dual problem is still separable with respect to both the units I 
and the time subdivisions p and the optimal commands ^ for time subdivision p, node n, 
unit t are immediately given by the following formulas: 

/ K(e)^ 
ui =0 if 7r„Q > Xn,p or 7r„C£ < Xn,p and ae < —-3 

, '^n,p ~ '''T,n ^max ^n,P else. 

4.1.2 Data and simulation protocol 

The data used for the simulations are inspired from real data. We suppose that each time step 
is divided in L time subdivisions also called hourly posts. Our generation strategy will be tested on 
a set of 456 scenarios. To each scenario is associated a reaHzation of the infiows in the reservoirs, 
the unavailability rates of the thermal units and of the demand at each time subdivision of the 
year. From these scenarios we build 3 different trees. Each tree corresponds to a vision more or 
less difficult of the evolution of the infiows and the demand on the coming year. We will call those 
trees Easy tree. Median tree and Difficult tree with evident interpretation of the predictions of the 
demands and the infiows on those trees. The scenario trees are trees of depth 364 days, with 5227 
nodes. At each node, we know the demand vector of the demands for all the posts p of this node, 
the infiows for all the hydro reservoirs and the programmed unavailability rates of all the hydro 
and thermal units. There are L = 3 hourly posts per day. We use the following generation units: 
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• Eleven thermal units. Every thermal unit i is described by its (unitary) generation cost, its 
maximal and minimal power, the number of thermal groups and the probability ai that a 
group works. 

• Two independent hydro plants. Each hydro plant is connected to a different reservoir. We 
know the maximal stock (in GWh) of each reservoir, the initial stock of each reservoir and the 
maximal power (in MW) of each plant. The maximal stock of the biggest reservoir is around 
30 times that of the other reservoir. This explains why we will essentially be interested in 
the evolution of the biggest reservoir stock on the year. 

• An EJP contract of 22 days with maximal available power: Pj = 2467 MW. 

Before presenting the results it remains to explain how the covariance matrix Q involved in VaRpA 
method is chosen. 



4.1.3 Calibration of the covariance matrix Q 

Since all we had were the scenarios of demands and the demands at the nodes of the three 
different trees generated from those scenarios following the Hues of 0, it was difficult to caHbrate 
the matrix Q. However, to have an idea of the impact of this method on the simulation process 
we supposed the demand at the different hourly posts were uncorrelated. We thus had to deem 
the diagonal elements of Q corresponding to the variances tT^(n,p) of the demand for node n and 
post p. This node n is associated to a time step t. To estimate a{n,p) we sort the demands of 
this time step by increasing order. Let's denote by d*--,, 1 <i < m the sample of demands ordered 
by increasing order for time step t. We then choose for a{i) which determines the uncertainty we 
have on d* .-j : 

a{i) = nim(^ '- ^ ' -), 

with the convention = and d(,„_|_i) = 2c?J„ — We will compare the results we obtained 

with the nominal model and with the VaR approaches VaRpA and VaRtenef- The outputs we 
are interested in are guided by the defaults of the nominal model outlined in subsection 12.21 We 
are thus interested in the distribution of the simulated costs and in the behavior of the hydro 
reservoirs. 



4.2 Numerical results 

4.2.1 Central and dispersion characteristics of the costs. 

We provide the mean and the standard deviation of the simulated costs. We also give the 
empirical quantile of order 0.95 (VaR 5%) and of order 0.99 (VaR 1%) of the distribution of these 
costs. We give the results using the three different trees (Easy, Median and Difficult tree) to solve 
the optimization problem and for all the methods. We also take a look at the method (called Mixt) 
consisting in cumulating the two modifications proposed by the two VaR approaches. 



Output 


Nominal 


Robust VaRpA 


Robust VaRbenef 


Mixt 


Mean 
St. Dev. 


467 529 291 
47 864 238 


488 271 561 
46 609 492 


459 515 733 
31 750 330 


459 991 109 
30 597 481 


VaR 1% 
VaR 5% 


671 599 214 
543 786 128 


672 311 733 
574 245 712 


557 798 070 
517 836 912 


558 719 609 
518 856 055 



Central and dispersion characteristics of the empirical distribution of the simulated 
costs on the Easy tree (management horizon of 1 year). 



11 



Output 


Nominal 


Robust VaRpA 


Robust VaRbenef 


Mixt 


Mean 
St. Dev 


466 498 435 
46 540 234 


486 587 235 
47 075 591 


462 184 762 
29 154 061 


462 334 138 
28 857 445 


VaR 1% 
VaR 5% 


689 406 053 
548 912 439 


679 536 989 
573 426 666 


557 838 702 
516 533 109 


554 297 141 
515 210 039 



Central and dispersion characteristics of the empirical distribution of the simulated 
costs on the Median tree (management horizon of 1 year). 



Output 


Nominal 


Robust VaRpA 


Robust VaRbenef 


Mixt 


Mean 
St. Dev. 


464 712 495 
44 615 708 


479 773 359 
47 172 746 


465 886 909 
28 146 821 


465 097 476 
28 272 953 


VaR 1% 
VaR 5% 


667 480 170 
543 728 341 


693 735 355 
562 040 194 


556 573 568 
517 110 690 


556 464 708 
518 142 699 



Central and dispersion characteristics of the empirical distribution of the simulated 
costs on the Difficult tree (management horizon of 1 year). 

We will essentially notice the following points: 

• The avcrago managing costs on the whole scenarios and for each of the three trees arc quite 
close for the four methods considered. For method VaRpA, the average costs are between 
3.2% and 4.4% greater than the nominal method. For method VaRbenef and Mixt, the costs 
can be greater or less than those of the nominal method. Nevertheless, those costs only vary 
between 0.08% and 1.7% compared with those of the nominal model. 

• Method VaRpA does not go into the good sense as the standard deviation of the costs as well 
as the VaR at 1% and 5% are greater than the same quantities computed for the nominal 
model. As for the methods VaRbenef and Mixt they lead in all the cases to reductions of 
the standard deviation of the costs (till 38% of reduction on the Difficult tree) and of the 
VaR at 1% and 5%. A reason that could explain the bad results of method VaRpA would 
be that the relaxation of the demand constraint in an ellipsoid works as an opportunity for 
the system to have an another reserve to perform its optimisation. But, this reserve does not 
exist in the Monte-Carlo simulation and thus, the strategy reveals itself to be too optimistic. 
The good results of VaRbenef could be attributed to the fact that the thermal problem is 
the optimization problem really dimensioning (the total thermal costs are around 100 times 
greater than the maximal valorization possible of the biggest reservoir) . It thus indeed seems 
interesting to envisage a robust approach which takes into account the only random involved 
in the thermal subproblem : the unavailability rates of the thermal plants. 



12 



4.2.2 Study of the trajectories of the biggest reservoir. 




Easy tree Median tree Difficult tree 

Figure 1: Evolution of the average stock (in MWh and on the whole scenarios) of the biggest 
reservoir during the year for all the methods and using the Easy, Median and Difficult trees as a 
support of the optimization process. 



For a given method, the strategy of management of the biggest reservoir softly vary when we 
change the scenario tree. The nominal and VaRpA methods tend to empty more the reservoir 
whose level increases at the end of the year. The methods VaRbenef and Mixt do not use the 
reservoir or very little at the beginning of the year. Globally, the reservoir has a higher level with 
those methods and is nearly full (its maximal stock is 3500 GWh) at the end of the year. In what 
follows, we say that the biggest reservoir (res for short) is at a low level if it contains at most 5% 
of its maximal stock. We say that a reservoir is at a high level if it attains a level greater or equal 
to its level at the beginning of the year less 5% of its maximal stock. Using those notations the 
array below permits to precise tendencies already observed on the above curves. 



# of weeks 




Hi 


gh level res 






Low level res 






Init 


VaRp 


A VaRbenef 


Mixt 


Init 


VaRpA 


VaRbenef 


Mixt 


1 


437 


405 


456 


456 


426 


456 


5 


9 


2 


421 


387 


456 


456 


423 


456 


4 


4 


3 


408 


377 


456 


456 


423 


456 


3 


3 


4 


390 


344 


456 


456 


417 


456 


3 


2 


5 


375 


312 


456 


456 


412 


456 


3 


2 


10 


201 


64 


456 


449 


365 


456 








15 


69 





449 


410 


256 


456 








20 


29 





438 


373 


92 


453 








25 


10 





414 


329 


14 


257 








30 


5 





371 


268 





25 









Number of scenarios among 456 for which the biggest reservoir is at least X weeks with a high or 
low level (optimizer launched on the Difficult tree and with all the methods). 
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4.2.3 Comparison of the distribution of the costs. 




Nominal Robust VaRpA 




3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 B.5 3.5 4 4.5 5 5.5 6 6.5 7 7.5 B 3.5 

Cost ^ ,q3 Cost ^ .njO 



Robust VaRbenef Robust Mixt 

Figure 2: Empirical densities of the management costs on the Difficult tree and for each method 
(management horizon 1 year). 



The empirical densities of the management costs for the nominal and VaRpA methods have 
tails of distribution bigger than those of methods VaRbenef and Mixt. A few scenarios are of very 
high cost for the nominal and VaRpA methods. On the contrary, the dispersion of the costs for 
the models VaRbenef and Mixt is smaller. We can illustrate those words by a few figures: 

• The scenario of highest cost corresponds to costs of 7.448 * 10^,8.232 * 10^,5.738 * 10^ and 
5.753 * 10* for respectfully the nominal VaRpA, VaRbenef and Mixt methods. 

• The less costly scenario for each of the models have close costs that are worth 3.822*10^,3.868* 
10^,3.997*10* and 3.991*10* for respectfully the nommal,VaRFA,VaRbenef and Mixt meth- 
ods. 

The shape below of the empirical cumulative distribution function of the management costs con- 
firms this tendency. The nominal method has a bigger density of scenarios whose cost is less than 
4.7* 10*. On the other hand, there exists a non neglectable number of scenarios of high costs. For 
VaRpA, it is worse as the density of scenarios whose cost is less than 4.7 * 10* is small. 
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Lower part Higher part 



Figure 3: Empirical cdf of the costs on the Difficult tree for the four methods. 



5 Conclusion 

This paper has presented an application of Value- At-Risk methods to robustify a stochastic 
optimization problem of yearly electric generation management. The starting point of our in- 
vestigation was a nominal model which conducted to a big standard deviation of the costs and 
which emptied too much the hydro reservoirs. Two models have been proposed to reduce those 
drawbacks. If the first model has not been concluding in practice, the second model has revealed 
very adapted to the objectives. On the one hand, this model conducts to a diminution of the 
standard deviation of the simulated management costs. On the other hand it tends to empty less 
the biggest reservoir. The success of this model comes from the fact that it modifies the thermal 
problem which is the optimizing problem really dimensioning. The model VaRpA, as for itself, 
is a classical robustification of the dual problem which takes into account the uncertainty on the 
demand. From a theoretical point of view, it permits to somehow stabilize the Lagrange multipliers 
which correspond to electricity prices in our application. On the other hand, this method should be 
interesting on difficult scenarios. To improve the results of this method, the problem of estimation 
of the matrix T should be carefully studied. 

Acknowledgement The authors are grateful to Anatoli luuditski of the Laboratoire de Mod- 
elisation et Controle de I'Universite Joseph Fourier for numerous advices and helpful discussions. 

A Appendix : Plant modeling and simulation process 

We first explain how the evolution of random variables is represented over the year. We then 
describe three important models of power generation units (thermal, hydro and demand side man- 
agement contracts called E JP) . We then briefiy comment the resolution of the different dual sub- 
problems. In subsection lA. 41 the way of determining a production schedule is outlined. 

A.l Prediction of random events and global problem 

Several stochastic optimization problems could be envisaged to model the problem of yearly 
electric generation management (e.g. |12[ll3[[n] V The possible evolutions of the realization of the 
random variables are represented by a markov chain (a set of scenarios organized in a tree) . Each 
node of the tree corresponds to a one day period of time. A day is divided in L hourly posts {Vn 
is the set of hourly posts of node n) . At each node n of the tree are attached a realization of the 
random variables we take into account : 
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• \n,p is the price of electricity for post p. 

• Vnip) is the electricity demand for post p. 

We will also use the following notations to describe the scenario tree: 

• T{n) is the time step associated to node n. 

• F{n) is the father node of node n. 

• S{n) is the set of son nodes of node n. 

• 7r„ is the probability to be at node n ( 7r„ = 1 V t) and TTT{n) is the probability to go 

T(n)=t 

from F{n) to n. 

• O is the set of nodes ( N = card O). 

• T is the last time step and Ot is the set of leaves. 

• dn^p is the duration of post p at node n. 

Other variables are attached to a node. They will be introduced in the description of the power 
generation units models. We give below the example of a scenario tree. In this example we have 
F(1) = 0, 5(2) = {4,5},... 




J Time(days) 




1 2 



Figure 4: Tree representing the different scenarios. 

A scenario in the tree is thus a path from the root node to a leaf node. The construction 
of the tree is based on aggregation procedures using historical data (demand, inflows,...). Those 
problems are discussed in |H].|rf5]. 

Given this representation of random events, the global problem of yearly power management 
scheduling consists in minimizing the average generation cost over the random tree while satis- 
fying the demand constraint and the operating constraints of the generation units. It expresses as 
Q . The reader should be aware that the solutions of this problem are indexed by the nodes of the 
scenario tree. If the scenario that occurs is represented in the tree we will have a generation sched- 
ule for this scenario. For a generation schedule that is not represented in the tree see subsection 

m 

A. 2 Modelling of power generation units 

Three kinds of generation units are modelled : the thermal, hydro and EJP units. 
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Thermal Units. Let i be the index of a thermal unit. The thermal units are completely de- 
scribed by two characteristics: 

1. The generation levels :-P^,p must remain between f^J™'" = and -P^'™^^ ~ '^'f n'^T n -'^max '^n,? 
with 

• „ the random unavailability rate, realization of a random variable r^. 

• t|,„ the programmed unavailability rate (deterministic). 

• -Pmax the maximal power (in MW) of unit t. 

2. The generation costs : O the thermal costs are 

where q is the unitary production cost for unit I. 

Hydro Units. An hydro valley is a set of interconnected plants and reservoirs with natural 
inflows in each reservoir. The state variables are the contents of the reservoir and the command 
variables are the discharge of the turbines and the water poured out of the reservoirs. In this yearly 
model, the total hydro production is aggregated in two different reservoirs non-interconnected. 
They both represent the total production capacity of the hydro generation units. The constraints 
are of two kinds: bounds on the volume of each reservoir, on the discharged water, and flow balance 
equations at each reservoir. Let I be an hydro plant. We will use the following notations: 

• is the content (in MWh) of the unique reservoir associated to i at the beginning of time 
step T{n). x^in and x^ax ^"^^ the lower and upper bounds (in MWh) on reservoir I level. 

• T^j ^ is the programmed unavailability rate (deterministic) 

• p is the natural inflow (in MWh) in reservoir t for node n, post p. = X]pep„ '^n,p- 

• -^max the maximal power of hydro plant £. 

• devf^ p is the amount of water (in MWh) poured outside of reservoir i and (;f, p is the 
discharge of plant £ (in MWh). vf^^™ = and vfi^'^^ are the lower and upper bounds 
on the discharged power. Thus the command variables for the hydro subproblems are 
«p)„,p = «p.^e<p)^_^. 

• Vff{x) is the value of the water stock x of reservoir i at the last time step T. The operating 
costs for all reservoirs are null for r(n) < T. Hence, only the value of the water stock at time 
T is taken into account. 

If jCh is the set of hydro units, the hydro problem thus consists in minimizing 
- ^ V^{xi + 4 - Z vi^p + devi^p), 



under the constraints 



^i=X%{n)+ H (4(n),p - 4(n),p - C«e«F(„),p) V {n,(.) 

xLin < 4 + X] ^°'n,p - vi,p - devi^p) < xi,^ y TIGOt 

pev„ (30) 

< vi^p < T^H,n PLx dn,p V {n,p, £) 

O^devi^^ ^ V {n,p,£) 

^min — — ^max ^ (^j-^)- 
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EJP contracts. An EJP contract I is represented by a production unit with the following 
features: 

• J I is the total number of days the contract can be used. 

• Each day, either the contract is used all day long or it is not used. A command variable 
defined for all node n permits to know whether the contract I is used at node n {tf^ = 1) or 
not (il =0). 

• is the stock (in days) still available on the contract £ for node n at the beginning of time 
step T{n). 

• Vj{.) is the function defining the value of an EJP stock at the last time step T for contract 

e. 

• The power linked to contract £ is Pj. 

Given a starting stock on EJP contract £ we maximize the value of the EJP stock at time step 
T which yields to the following problem: 

ti) 

for the root node. 

VtigO (31) 
V (n, £) 
y{n,£) 
\fneOT, V£, 

if £j is the set of EJP contracts. Thus if ^ S £j, the EJP command ^ is given by i„ Pj dn^p- To 
be complete, we should take into account nuclear power plants. However, this would not change 
much things as they can be modelled in a similar way (see [El)- Using this modelling for the 
generation units, the constraint of satisfaction of the demand writes Vnijp) — ^ Ai,p(^n,pi ^l,p) = 



si>0 
si~ti>0 



A. 3 Solving the different dual subproblems 

We briefiy detail the first step of the space decomposition algorithm given in subsection 14.1.11 
The nominal thermal dual subproblem has an evident solution and the hydro subproblem is a 
linear optimization problem of big size (around 36 000 variables in our case) solved using interior 
point methods. The EJP problem is an NP complete, non convex optimization problem solved 
using dynamic stochastic programming. We know, for every contract £, the values at the last time 
step T of all possible values of the EJP stock x : Vj{x,T). We deduce, using HJB equations, 
(backward phase) the bellman values V^{x,n) for all contract £ and all node n: 

^ ' ^"l x-ti>0, tie {0,1}. ^^'^ 

Knowing the stock of every contract £ at the beginning of the year, the forward phase consists in 
deducing the optimal EJP commands to apply using the bellman values computed in the backward 
phase. Both kinds of methods used (dynamic stochastic programming and interior points methods) 
have a complexity that depends on the dimension of the dual space . A known drawback of dynamic 
programming is that its complexity grows exponentially with the state variable dimension. 
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A. 4 The simulation process 



The resolution of (P) provides optimal marginal prices A* ^ that are useful to elaborate a 
strategy which, for any realization of the random variables on the time period, will allow the 
computation of a generation schedule. This strategy has the form of Bellman functions and are 
computed with the fohowing version of the Bellman principle. Let a reserve i be given. For the 
last time step T, the Bellman function V^{x, n) is known for each stock level x and node n of the 
leaves. Between two grid points, the value is supposed to be hnear. The Bellman values V^{x,n) 
for all the nodes n of the tree and all stock step x are computed using the following recursive 
formula: 

y^(a;,n) = max( ^ 7rT(m)(y^(a; + - ^ vi^p + devi^p,m) 
mes{n) per„ 

x + ai- xi,^^ < vi^p + devi^p <x + ai 
peVn 

[ < vi^p < < devi^p 



(33) 



if £ stands for an hydro reservoir and 



y^(x,n)=max( ^ M^Wix - Cm) + dn,pPXA)) 

meS(n) pev„ 
ti{ti-l) = 0, x-ti>0 



(34) 



if £ stands for an EJP contract. The Bellman function for time step t is then given by: V^{x,t) = 
Sne t ^{njV^lx,!!,) for stock x, unit i. 

This algorithm is only a stochastic dynamic programming (SDP) performed on marginal values. 
Once those Bellman functions are computed, it is possible to perform a Monte-Carlo simulation of 
the generation scheduling using 5x^^ (^x^t) as a "fuel cost" of the energy kept in the reserve £. 
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